program case_selector

    use iasing_osse
    use netcdf

    implicit none
    ! condition for optional parameters
    logical :: lopt_metview = .true.

    ! Counters
    integer :: jo

    ! Arguments
    integer :: nargs
    character(len=64) :: cargs(10)

    ! Netcdf vars
    character(len=64) :: dbfile
    integer :: dbid, natm, nlev
    type(r4a_ids) :: db_ids

    ! Netcdf arrays for criteria selection
    real(JPR4), dimension(:), allocatable    :: lat, lon, lsm, sza, ele
    real(JPR4), dimension(:), allocatable    :: iaa, iza, saa, tsu, psu
    integer(JPI4), dimension(:), allocatable :: clc, dat, tim
    logical, dimension(:), allocatable :: lcond
    integer, dimension(:), allocatable :: jsel, indexes
    integer :: nsel

    ! Includes from libiasing_osse
    include "r60_initializer.interface"
    include "checknc.interface"

    ! Get arguments
    nargs = command_argument_count()
    if (nargs.eq.1) then
        call get_command_argument(1,dbfile)
    else
        print*, "No file provided or bad number of arguments. STOP"
        stop -1
        !dbfile = "/cnrm/obs/data2/andreyj/4ABis_runs/db_20130806_noise.nc"
    endif

    ! Open and initilise databank
    call checknc(nf90_open(dbfile, NF90_NOWRITE, dbid))
    db_ids = r60_initializer(dbid)

    ! Inquire number of atmospheres
    call checknc(nf90_inquire_dimension(dbid, db_ids%atm_did, len=natm))
    call checknc(nf90_inquire_dimension(dbid, db_ids%lev_did, len=nlev))
    print*, "dbfile: ", trim(dbfile)
    print*, "natm:   ", natm
    print*, "nlev:   ", nlev
    print*

    ! Allocate arrays for the different parameters
    ! Coordinates
    allocate(lon(natm))
    call checknc(nf90_get_var(db_ids%coo%gid, db_ids%coo%lon, lon))
    allocate(lat(natm))
    call checknc(nf90_get_var(db_ids%coo%gid, db_ids%coo%lat, lat))
    allocate(ele(natm))
    call checknc(nf90_get_var(db_ids%coo%gid, db_ids%coo%ele, ele))
    allocate(dat(natm))
    call checknc(nf90_get_var(db_ids%coo%gid, db_ids%coo%dat, dat))
    allocate(tim(natm))
    call checknc(nf90_get_var(db_ids%coo%gid, db_ids%coo%tim, tim))
    ! Observation parameters
    allocate(psu(natm))
    call checknc(nf90_get_var(db_ids%opa%gid, db_ids%opa%psu, psu))
    allocate(tsu(natm))
    call checknc(nf90_get_var(db_ids%opa%gid, db_ids%opa%tsu, tsu))
    allocate(iza(natm))
    call checknc(nf90_get_var(db_ids%opa%gid, db_ids%opa%iza, iza))
    allocate(iaa(natm))
    call checknc(nf90_get_var(db_ids%opa%gid, db_ids%opa%iaa, iaa))
    allocate(sza(natm))
    call checknc(nf90_get_var(db_ids%opa%gid, db_ids%opa%sza, sza))
    allocate(saa(natm))
    call checknc(nf90_get_var(db_ids%opa%gid, db_ids%opa%saa, saa))
    allocate(lsm(natm))
    call checknc(nf90_get_var(db_ids%opa%gid, db_ids%opa%lsm, lsm))
    allocate(clc(natm))
    call checknc(nf90_get_var(db_ids%opa%gid, db_ids%opa%clc, clc))

    ! Arrays for logical condition
    allocate(lcond(natm))
    allocate(jsel(natm))
    allocate(indexes(natm))
    indexes=(/(jo, jo=1,natm)/)

    ! Read array values

    ! Close netcdf file
    call checknc(nf90_close(dbid))

    ! Apply condition ...
    lcond = lsm.eq.0 .and.  &
        clc .eq. 0   .and.  &
        lon .gt. -60   .and.  &
        lon .lt. -15   .and.  &
        lat .gt.30.  .and.  &
        lat .lt.60. 
    !lcond = lsm == 0 .and.  &
    !    clc == 0   .and.    &
    !    abs(lat) >= 70.
    !lcond = .true.
    !  ... and count number of cases
    nsel = count(lcond)
    jsel(1:nsel) = pack(indexes,lcond)

    ! Print information
    print*, "nsel:   ", nsel
    print*

    ! Write results
    open(10, file="selected")
    write(10,'(A)') dbfile
    write(10,'(I8)') nsel, jsel(1:nsel)
    close(10)

    ! Write geopoints file for metview
    if (lopt_metview .eqv. .true.) then
        open(11, file="selected.metview", action='write')
        write(11,'(2(A,/),A)') "#GEO", "#FORMAT LLV", "#DATA"
        do jo = 1, nsel
            write(11,'(2F10.3, I8)') lon(jsel(jo)), lat(jsel(jo)), 0
        enddo
        close(11)
    endif
    
    ! Deallocate arrays
    deallocate(lon)
    deallocate(lat)
    deallocate(lsm)
    deallocate(sza)
    deallocate(clc)
    deallocate(lcond)
    deallocate(jsel)
    deallocate(indexes)


end program
